{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a9cac518-0f99-4df6-9d13-0a0b8aa54827",
   "metadata": {},
   "source": [
    "# Dirichlet problem in a rectangular domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b5d9a45-390c-4cb1-965e-89ee51b1a635",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.special as sc\n",
    "import calfem.geometry as cfg\n",
    "import calfem.vis as cfv\n",
    "\n",
    "def create_geometry(el_size=0.2):\n",
    "    g = cfg.Geometry()\n",
    "\n",
    "    points = [[0.0, 0.0],[1.0, 0.0],[1.0, 1.0],[0.0, 1.0],[0.5,0.5]]\n",
    "    splines = [[0, 1],[1, 2],[2, 3],[3, 0]]\n",
    "    surfaces = [[0,1,2,3]]\n",
    "    \n",
    "    for x, y in points:\n",
    "        g.point([x, y], el_size=el_size)\n",
    "\n",
    "    for s in splines:\n",
    "        g.spline(s)\n",
    "\n",
    "    for f in surfaces:\n",
    "        g.surf(f)\n",
    "        # g.struct_surf(f)\n",
    "\n",
    "\n",
    "    bcs = {101:[0],\n",
    "           201:[1],\n",
    "           301:[2],\n",
    "           401:[3]\n",
    "           }    \n",
    "    \n",
    "    for marker in bcs:\n",
    "        for i in bcs[marker]:\n",
    "            g.curve_marker(ID=i, marker=marker)    \n",
    "    \n",
    "    # cfv.drawGeometry(g)\n",
    "    # cfv.showAndWait()\n",
    "    \n",
    "    return g\n",
    "\n",
    "def neu101(x,y):\n",
    "    return -exact_uy(x,y)\n",
    "\n",
    "def neu201(x,y):\n",
    "    return exact_ux(x,y)\n",
    "\n",
    "def neu301(x,y):\n",
    "    return exact_uy(x,y)\n",
    "\n",
    "def neu401(x,y):\n",
    "    return -exact_ux(x,y)\n",
    "\n",
    "########################## paper example 2a #################################\n",
    "# def source(x,y):\n",
    "#     return np.ones_like(x)\n",
    "\n",
    "# def exact_u(x,y):\n",
    "#     return x*y-3*x*y**2-x**2/2+x**3\n",
    "\n",
    "# def exact_ux(x,y):\n",
    "#     return y-3*y**2-x+3*x**2\n",
    "\n",
    "# def exact_uy(x,y):\n",
    "#     return x-6*x*y\n",
    "\n",
    "# def q0_x(x,y):\n",
    "#     return -x\n",
    "\n",
    "# def q0_y(x,y):\n",
    "#     return np.zeros_like(x)\n",
    "\n",
    "#################################################################\n",
    "# def source(x,y):\n",
    "#     return 2*(x-x**2)+2*(y-y**2)\n",
    "\n",
    "# def exact_u(x,y):\n",
    "#     return x*(1-x)*y*(1-y)\n",
    "\n",
    "# def exact_ux(x,y):\n",
    "#     return (1-2*x)*y*(1-y)\n",
    "\n",
    "# def exact_uy(x,y):\n",
    "#     return x*(1-x)*(1-2*y)\n",
    "\n",
    "# def q0_x(x,y):\n",
    "#     return -(2*(x**2/2-x**3/3)+2*(y-y**2)*x)\n",
    "\n",
    "# def q0_y(x,y):\n",
    "#     return np.zeros_like(x)\n",
    "############################## paper example 2b ###################################\n",
    "# xc = -0.5\n",
    "# yc = -0.5\n",
    "# def source(x,y):\n",
    "#     return -1/np.sqrt((x-xc)**2+(y-yc)**2)\n",
    "\n",
    "# def exact_u(x,y):\n",
    "#     return np.sqrt((x-xc)**2+(y-yc)**2)\n",
    "\n",
    "# def exact_ux(x,y):\n",
    "#     return (x-xc)/np.sqrt((x-xc)**2+(y-yc)**2)\n",
    "\n",
    "# def exact_uy(x,y):\n",
    "#     return (y-yc)/np.sqrt((x-xc)**2+(y-yc)**2)\n",
    "\n",
    "# def q0_x(x,y):\n",
    "#     return np.log(x-xc+np.sqrt((x-xc)**2+y**2-2*y*yc+yc**2))\n",
    "\n",
    "# def q0_y(x,y):\n",
    "#     return np.zeros_like(x)\n",
    "############################## reference paper  ############################\n",
    "\n",
    "def source(x,y):\n",
    "    return 2*np.pi*np.cos(2*np.pi*y)*(2*np.pi*(1+x**2)*np.cos(np.pi*x**2)+np.sin(np.pi*x**2))\n",
    "\n",
    "def exact_u(x,y):\n",
    "    return np.cos(np.pi*x**2)*np.cos(2*np.pi*y)\n",
    "\n",
    "def exact_ux(x,y):\n",
    "    return -2*np.pi*x*np.sin(np.pi*x**2)*np.cos(2*np.pi*y)\n",
    "\n",
    "def exact_uy(x,y):\n",
    "    return -2*np.pi*np.cos(np.pi*x**2)*np.sin(2*np.pi*y)\n",
    "\n",
    "def q0_x(x,y):\n",
    "    return -(2*np.sqrt(2)*np.pi**2*np.cos(2*np.pi*y)*sc.fresnel(np.sqrt(2)*x)[1]+2*np.pi*x*np.cos(2*np.pi*y)*np.sin(np.pi*x**2))\n",
    "\n",
    "def q0_y(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "\n",
    "\n",
    "import toolkits\n",
    "from importlib import reload\n",
    "from fem_max import Poisson2D\n",
    "\n",
    "geometry = create_geometry(el_size=0.1)\n",
    "Dirichlet = []\n",
    "Neumann = [101,201,301,401]\n",
    "hole = []\n",
    "BCfunc = {101:neu101, 201:neu201, 301:neu301,401:neu401}\n",
    "\n",
    "model = Poisson2D(geometry, Dirichlet, Neumann, hole, BCfunc, source, q0_x, q0_y)\n",
    "model.mesh(display=False)\n",
    "# model.efem(display=False)\n",
    "# model.fem(display=False)\n",
    "# model._display_efem(cmap=\"jet\")\n",
    "# model._display_fem(cmap=\"jet\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EFEM",
   "language": "python",
   "name": "efem"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
